Truncated-Determinant Diagrammatic Monte Carlo for Fermions with Contact 
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For some models of interacting fermions the known solution to the notorious sign-problem in Monte 
Carlo (MC) simulations is to work with macroscopic fermionic determinants; the price, however, 
is a macroscopic scaling of the numerical effort spent on elementary local updates. We find that 
the ratio of two macroscopic determinants can be found with any desired accuracy by considering 
truncated (local in space and time) matices. In this respect, MC for interacting fermionic systems 
becomes similar to that for the sign-problem-free bosonic systems with system-size independent 
update cost. We demonstrate the utility of the truncated-determinant method by simulating the 
attractive Ifubbard model within the MC scheme based on partially summed Feynman diagrams. We 
conjecture that similar approach may be useful in other implementations of the sign-free determinant 
schemes. 
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Monte Carlo (MC) methods are a unique tool for 
studying large interacting systems. The most severe lim- 
itation on their applicability is imposed by the so-called 
sign-problem (SP) when relevant contributions to statis- 
tics alternate in sign and almost exactly compensate each 
other in the final answer . Frustrating interactions and 
anticommutation relations for fermion operators are typ- 
ically at the origin of the sign-problem. In this Letter, we 
address the case when quantum statistics is non-positive 
only because of the fermion exchange cycles. 

One solution to the fermion SP is offered by the de- 
terminant Monte Carlo (DetMC) (see, e.g. The 
idea is that all contributions to the many-body statis- 
tics obtained by exchanging fermion places (in a cer- 
tain representation) can be written as a product of two 
determinants — for spin-up and spin-down species — and 
in cases when the two real determinants coincide the re- 
sult is positive definite. In Metropolis-type algorithms 
0, MC updates are accepted with probabilities pro- 
portional to the ratio of final and initial configuration 
weights; in DetMC the corresponding acceptance ratio, 
i?, is based on the ratio of large determinants. Unfortu- 
nately, calculating determinants ratio for macroscopically 
large matrices is very expensive numerically: even with 
tricks involving the Hubbard-Stratonovich tranformation 
the algorithm proposed by Blankenbecler, Scalapino and 
Sugar 0, 0| still requires L^'^ operations per update for 
a c?-dimensional system with L lattice points per dimen- 
sion. The same scaling is true for the continuous-time 
scheme 6]. In contrast, for bosonic systems with local 
interactions the number of operations per update is small 
and system-size independent, i.e. they can be simulated 
V^"^ times faster! 

Since the bottleneck of DetMC is the calculation of 
i?, one may question the paradigm of calculating it "ex- 
actly" : In any case, computer operations always involve 
systematic round-off errors. The other example is pro- 



vided by (pseudo)random number generators — they are 
always imperfect and result in systematic errors equiva- 
lent to small errors in R which, however, remain practi- 
cally undetectable (for good generators) in final results. 
The heuristic explanation of why small errors in R do 
not ruin the simulation is as follows. The Metropolis al- 
gorithm is a scheme with strong relaxation towards equi- 
librium distribution, and local configuration updates may 
be viewed as the result of dissipative coupling to the ther- 
mal bath (this picture is often used to model dissipative 
kinetics |l|). Uncontrolled errors in R may then be re- 
garded as a small stochastic noise in the relaxational dy- 
namics. As such, it only slightly modifies the equilibrium 
state and its properties. This is a standard argument in 
the linear response theory. 

It seems natural then to suggest that if the goal is to 
simulate the result with n-digit accuracy, there is no need 
to calculate the acceptance ratio with accuracy much 
higher than n digits. Often, we simply ignore this issue 
because getting R with machine precision does not cost 
any extra CPU time. In determinant methods, however, 
there is a potential of huge efficiency gains if approximate 
values of R can be calculated much faster. We demon- 
strate the feasibility of this approach by showing that the 
ratio of two macroscopic determinants can be found with 
high accuracy by considering truncated matrices dealing 
only with the local (in space-imaginary time) structure 
of the configuration space. The computational cost of 
updates in the corresponding "truncated-determinant" 
scheme is system-size independent — an efficiency increase 
c>c L?'^ for large L. 

In what follows, we discuss the solution of the Hubbard 
model for fermions within a simple diagrammatic MC 
scheme based on Feynman diagrams partially summed 
over fermion propagator permutations 0, 0|- The re- 
sulting diagram weight is the square of the determinant 
composed of finite-temperature fermion propagators. We 
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explain how the determinant ratio for local updates may 
be calculated using truncated matrices, and demonstrate 
the feasibility of the proposed approach. We also show 
that going to larger system sizes has little effect on the 
scheme performance. Finally, we conjecture that large 
efficiency gains are expected in other sign-problem free 
DetMC schemes, e.g. in lattice QCD simulations with 
quark fields 0. It is also worth noting that in the dia- 
grammatic DetMC scheme the update cost does not de- 
pend directly on the lattice period, which is a big advan- 
tage for simulations of dilute systems. By "dilute" we 
mean dilute with respect to the lattice, but not neces- 
sarily with respect to the interaction; the latter can be 
effectively large due to a resonance on a (quasi)-bound 
state. Correspondingly, the new scheme is very promis- 
ing for the study of ultra cold fcrmionic systems in the 
regime of strong Feshbach resonant interaction, including 
the crossover from the Bardeen-Cooper-SchriefFer pairing 
to the Bose-Einstein condensation of molecules [i^- 

Model and method. We consider interacting lattice 
fermions with the Hubbard Hamiltonian H = Hq + Hi : 



Hn = 



ko- 



(eic - Ai)cLcka, Hi = C/^nxT^xi, (1) 



where c^.^ is the fermion creation operator, rtxo- = c^o-Cxct, 
cr = t, J, is the spin index, x runs over the L"* points of the 
simple cubic lattice, k runs over the corresponding Bril- 
louin zone, ek = — 2t ^^^^ cos fcafl is the tight-binding 
dispersion law, and /i is the chemical potential. For defi- 
niteness and numerical tests, we confine ourselves to the 
d = 2 spacial lattice with periodic boundary conditions. 
We use the hopping amplitude, t, and lattice constant a 
as units of energy and distance, respectively. 

Following Refs. 0, Q we start with writing the statis- 
tical operator in the interaction representation. 
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where Hi{t) — e^°'^Hie ^"'^ and Tr is the time ordering 
operator, and expanding it in powers of Hi: 



n „ Jt 



p=0 



X Tr 
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(3) 



This expansion for the partition function generates stan- 
dard Feynman diagrams jsl- Graphically, each term is a 
set of four-point vertices with two incoming (spin-up and 
spin-down), and two outgoing (spin-up and spin-down) 
lines which connect vertices. Each line is associated with 
the imaginary time fermion propagator, G(j(xi — Xj, — 
TJ;^i, P)^- Tr [r,e-'3«°c„(x,T,)4 (x,r,)] . 

A straightforward MC sampling of diagrammatic se- 
ries would be impossible because of the sign-problem. 



However, if for a given configuration of p vertices, Sp = 
{(xj,T,), j = 1, . . . one sums over all (p!)^ ways of 
connecting them by propagators, then the result can be 
written as a product of two determinants, one for spin up, 
and another for spin down (see e.g. 0)- The differential 
weight of the vertex configuration (or vertex diagram) is 
then 

p 

dV{Sp) = {-Uf detAT(5p)detA^(5p) J|dr, , (4) 



i=l 



where A'^{Sp) are p x p matrices: A° 



Gcr(Xi 



Xj, Ti — Tj). For equal number of up- and down-particles, 
det A^ det A^- = [det A]^, and negative U the vertex dia- 
gram weight is always positive. [At half filling, +ni = 
1, the sign of U changes when hole representation is used 
for one of the spin components, so this scheme may be 
also used for the repulsive Hubbard model.] 

MC sampling in the vertex configuration space ip,Sp) 
can be per formed by standard MC rules (see, e.g. 
Refs. [illlj) using just one pair of complementary up- 
dates "D and C: inV one selects at random one of the ver- 
tices and suggests to delete it from the configuration; in C 
an additional vertex is suggested to be inserted at some 
point randomly selected in the space-time box (] x L'^. 
These updates decrease/increase the rank of A by one. 
The acceptance ratio for the V/C pair of updates is then 
based on the ratio of two determinants 

Rp = dot A(5p+i)/ det A{Sp) , (5) 

where Sp+i = {Sp, (xp+i,Tp+i)} (we omit the spin index 
for brevity). 

The bottleneck of this simple scheme is in evaluating 
Rp when p is macroscopically large. The typical number 
of vertices is determined by the number of particles, in- 
teraction strength, and inverse temperature asp oc NpU. 
The truncated-determinant idea is to calculate Rp much 
faster at the expense of accuracy using the following con- 
jecture originating from physical, rather than mathemat- 
ical, arguments. The vertex configuration represents a 
sequence of virtual particle collisions in the many-body 
system, and it is likely that local changes in its structure 
depend only on the immediate neighborhood of the up- 
dated region. [We note that the idea of employing the 
local nature of the fermion-boson coupling has been used 
in but it has not been extended to the fermionic de- 
terminant.] Quantitatively, we define a norm, || . . . ||, or a 
distance, between vertices in space-time (several choices 
are discussed below), and construct a truncated vertex 
configuration, Sp^"^ , such that all points in Sp^"^ satisfy 



(Xp+l,Tp+l)|| < 



(6) 



Correspondingly, S^'^li = {5^^', (xp+iTp+i)}. We may 
now use truncated configurations to calculate the ratio 
(ISJ approximately as 



i?W=detA(4^^^i)/detA(4^)). 



(7) 
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Clearly, when £ ^ L we recover the exact ratio. Our 

(£) 

conjecture is then that Rp quickly converges to Rp and 
there exists a healing length in the (x, T)-space charac- 
terizing this convergence. If this is the case, then £ may 
be considered as a microscopic (system-size independent) 
parameter controlling the accuracy and cfhciency of sim- 
ulation. 

The proper choice of the norm || . . . || depends on sys- 
tem parameters. In the strongly correlated case the nat- 
ural units of distance and time are provided by the Fermi 
momentum, kp, and Fermi energy ep- One possibility is 
then 



(x,r)-(x',r') 



A:|,(x-x')2+e|,(r-r')2 . (8) 



Geometrically, this measure results in a set of ver- 
tices Sp inside the space-time ellipsoid centered at 
(xp+i,rp+i). For dense systems ||(x,t) - (x',r')|| = 
max{|x — x'l, |t — r'|} is equally appropriate and our 
data for the largest system were obtained using this mea- 
sure. At temperatures comparable to ep one may ac- 
count for all vertices in the f-direction and simply write 
||(x, r) — (x', r')||cyi = kp\K — The corresponding 
geometrical figure is a /3-cylinder. Similarly, in small sys- 
tems one may consider truncating configurations only in 
time direction. 

Numerical results and discussion. Our tests of the 
truncated-determinant scheme were done for the attrac- 
tive Hubbard model with U = —4, ^ = —2, /3 = 10, 
and periodic boundary conditions. First, we simulated 
a small = 4^ cluster for which the ground state en- 
ergy (of 10 particles) is known from exact diagonalization 
studies 01 . Since the spatial dimension is so small we 
truncate vertex configurations only in the imaginary time 
direction. In Fig. ^ we show how the result for energy 
converges to the exact value. We stress, that at all stages 
of the MC simulation we never even write the full config- 
uration determinant, which rank is about 2.5 times larger 
than typical values of p for the cutofi^ radius 2. 

In Fig. 121 we present our data for the = 100^ 
system — now, using determinant truncation both in 
space and in time directions, Eq. 0. Remarkably, the 
convergence is achieved around the same value of the 
truncation radius, which proves that the computational 
cost per update is not subject to macroscopic scaling. 
It is instructive to see how data convergence for energy 
correlates with the typical errors introduced by the ap- 
proximate calculation of Rp. In Fig. we show examples 
of Rp dependence on the truncation radius for a number 
of randomly selected MC configurations. Clearly, quite 
large fluctuations in Rp are statistically "averaged out" 
in the final result for energy. We are not aware of any 
other method capable of simulating fermionic systems of 
comparable size. 

Apart from the Hubbard model tests, we have also ver- 
ified that the use of truncated-determinants for randomly 
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FIG. 1: Energy and density dependence on the imaginary 
time cutoff for the Hubbard model for L — 4. Dots represent 
the MC data, and solid lines are the exact diagonalization 
results for 10 particles [TH . 
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FIG. 2: Potential energy dependence on the truncation radius 
for the Hubbard model for L — 100. 



seeded vertex configurations works as nicely to speed up 
the calculation of Rp. 

The benchmark DetMC method by Blankenbecler, 
Scalapino and Sugar (BSS) 0, 01 is based on the 
Trotter-Suzuki imaginary time slicing, and the Hubbard- 
Stratonovich transformation 0| . The rank of the matrix 
used in the calculation of the acceptance ratio equals the 
number of lattice sites L'', and ~ L^*^ operations are re- 
quired to find its determinant. The necessity of handling 
large matrices, although made possible by this method, 
requires both elaborate finite-size scaling analysis of MC 
data 01, and special efforts for the calculation stabi- 
lization at low temperatures The contour-distortion 
stabilization techniques (see, e.g. [l^ fisj') help to alle- 
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FIG. 3: Determinant ratios Rp as functions of £ for randomly 
selected MC configurations for the L = 100 system. 



viate the sign problem, but still suffer from the severe 
scaling of the computational cost per update. More re- 
cently, Rombouts, Heide, and Jachowicz |y] improved the 
BSS scheme by formulating it in the r-continuum. ft is 
easy to directly compare this scheme with ours because 
the starting point is exactly the same — the expansion of 
the statistical operator in powers of U. Rombouts et al. 
used the auxiliary Ising variables to decompose four-point 
vertices into the sums of single-particle exponentials, and 
arrived at the number of operations for performing one 
update scaling as L^'^. At this point we notice that while 
we work with the same vertex configuration structure, in 
our scheme there is no extra summation over the auxil- 
iary variables, and the calculation of the acceptance ratio 
is system-size independent. 

Recently, a substantial improvement in lattice QCD 
simulations has been achieved by including quark-loop 
effects (see, e.g. 0). After the quarks are integrated 
out, their effects are described by the macroscopic posi- 
tive definite determinant, det A([/), where U is the con- 
figuration of gluon matrices. The conjecture is that the 
ratio det A(C/)/ det A(t/') for local updates of gluon ma- 
trices, U — > C/', can be calculated by accounting only for 
the immediate neighborhood of the updated lattice bond. 



We doubt that the truncated-determinant schemes will 
help to speed up simulations with the sign-problem. If 
the average configuration sign is small, the answer is de- 
termined by small differences between the sign-positive 
and sign-negative contributions, i.e. each contribution 
has to be calculated to much higher accuracy than would 
be sufficient in the positive definite case. As a result, 
higher and higher precision is required for Rp, and ad- 
vantages of our approach quickly vanish. 
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